// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1_CUBE3D_LOCALBASIS_HH
#define DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1_CUBE3D_LOCALBASIS_HH

#include <numeric>
#include <vector>

#include <dune/common/fmatrix.hh>

#include "../../common/localbasis.hh"

namespace Dune
{
  /**
   * \ingroup LocalBasisImplementation
   * \brief First order Raviart-Thomas shape functions on the reference hexahedron.
   *
   * \tparam D Type to represent the field in the domain.
   * \tparam R Type to represent the field in the range.
   *
   * \nosubgrouping
   */
  template<class D, class R>
  class RT1Cube3DLocalBasis
  {

  public:
    typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,3,Dune::FieldVector<R,3>,
        Dune::FieldMatrix<R,3,3> > Traits;

    /**
     * \brief Make set number s, where 0 <= s < 64
     *
     * \param s Edge orientation indicator
     */
    RT1Cube3DLocalBasis (unsigned int s = 0)
    {
      sign0 = sign1 = sign2 = sign3 = sign4 = sign5 = 1.0;
      if (s & 1)
      {
        sign0 = -1.0;
      }
      if (s & 2)
      {
        sign1 = -1.0;
      }
      if (s & 4)
      {
        sign2 = -1.0;
      }
      if (s & 8)
      {
        sign3 = -1.0;
      }
      if (s & 16)
      {
        sign4 = -1.0;
      }
      if (s & 32)
      {
        sign5 = -1.0;
      }
    }

    //! \brief number of shape functions
    unsigned int size () const
    {
      return 36;
    }

    /**
     * \brief Evaluate all shape functions
     *
     * \param in Position
     * \param out return value
     */
    inline void evaluateFunction (const typename Traits::DomainType& in,
                                  std::vector<typename Traits::RangeType>& out) const
    {
      out.resize(36);

      out[0][0] = sign0*(-3.0*in[0]*in[0] + 4.0*in[0] - 1.0);
      out[0][1] = 0.0;
      out[0][2] = 0.0;

      out[1][0] = sign1*(-2.0*in[0] + 3.0*in[0]*in[0]);
      out[1][1] = 0.0;
      out[1][2] = 0.0;

      out[2][0] = 0.0;
      out[2][1] = sign2*(-3.0*in[1]*in[1] + 4.0*in[1] - 1.0);
      out[2][2] = 0.0;

      out[3][0] = 0.0;
      out[3][1] = sign3*(3.0*in[1]*in[1] - 2.0*in[1]);
      out[3][2] = 0.0;

      out[4][0] = 0.0;
      out[4][1] = 0.0;
      out[4][2] = sign4*(-3.0*in[2]*in[2] + 4.0*in[2] - 1.0);

      out[5][0] = 0.0;
      out[5][1] = 0.0;
      out[5][2] = sign5*(3.0*in[2]*in[2] - 2.0*in[2]);

      out[6][0]  = -18.0*in[0]*in[0]*in[1] + 9.0*in[0]*in[0] + 24.0*in[0]*in[1] - 6.0*in[1] -12.0*in[0] + 3.0;
      out[6][1]  = 0.0;
      out[6][2]  = 0.0;

      out[7][0]  = -18.0*in[0]*in[0]*in[1] + 9.0*in[0]*in[0] + 12.0*in[0]*in[1] - 6.0*in[0];
      out[7][1]  = 0.0;
      out[7][2]  = 0.0;

      out[8][0]  = 0.0;
      out[8][1]  = 18.0*in[0]*in[1]*in[1] - 9.0*in[1]*in[1] - 24.0*in[0]*in[1] + 12.0*in[1] + 6.0*in[0] - 3.0;
      out[8][2]  = 0.0;

      out[9][0]  = 0.0;
      out[9][1]  = 6.0*in[1] - 12.0*in[0]*in[1] - 9.0*in[1]*in[1] + 18.0*in[0]*in[1]*in[1];
      out[9][2]  = 0.0;

      out[10][0] = 0.0;
      out[10][1] = 0.0;
      out[10][2] = -3.0 + 6.0*in[0] + 12.0*in[2] - 24.0*in[0]*in[2] - 9.0*in[2]*in[2] + 18.0*in[0]*in[2]*in[2];

      out[11][0] = 0.0;
      out[11][1] = 0.0;
      out[11][2] = 6.0*in[2] - 12.0*in[0]*in[2] - 9.0*in[2]*in[2] + 18.0*in[0]*in[2]*in[2];

      out[12][0] = 3.0 - 12.0*in[0] - 6.0*in[2] + 24.0*in[0]*in[2] + 9.0*in[0]*in[0] - 18.0*in[0]*in[0]*in[2];
      out[12][1] = 0.0;
      out[12][2] = 0.0;

      out[13][0] = -6.0*in[0] + 12.0*in[0]*in[2] + 9.0*in[0]*in[0] - 18.0*in[0]*in[0]*in[2];
      out[13][1] = 0.0;
      out[13][2] = 0.0;

      out[14][0] = 0.0;
      out[14][1] = 3.0 - 12.0*in[1] - 6.0*in[2] + 24.0*in[1]*in[2] + 9.0*in[1]*in[1] - 18.0*in[1]*in[1]*in[2];
      out[14][2] = 0.0;

      out[15][0] = 0.0;
      out[15][1] = -6.0*in[1] + 12.0*in[1]*in[2] + 9.0*in[1]*in[1] - 18.0*in[1]*in[1]*in[2];
      out[15][2] = 0.0;

      out[16][0] = 0.0;
      out[16][1] = 0.0;
      out[16][2] = -3.0 + 6.0*in[1] + 12.0*in[2] - 24.0*in[1]*in[2] - 9.0*in[2]*in[2] + 18.0*in[1]*in[2]*in[2];

      out[17][0] = 0.0;
      out[17][1] = 0.0;
      out[17][2] = 6.0*in[2] - 12.0*in[1]*in[2] - 9.0*in[2]*in[2] + 18.0*in[1]*in[2]*in[2];

      out[18][0] = -9.0 + 36.0*in[0] + 18.0*in[1] + 18.0*in[2] - 72.0*in[0]*in[1] - 72.0*in[0]*in[2] - 36.0*in[1]*in[2] + 144.0*in[0]*in[1]*in[2] - 27.0*in[0]*in[0] + 54.0*in[0]*in[0]*in[1] + 54.0*in[0]*in[0]*in[2] - 108.0*in[0]*in[0]*in[1]*in[2];
      out[18][1] = 0.0;
      out[18][2] = 0.0;

      out[19][0] = 18.0*in[0] - 36.0*in[0]*in[1] - 36.0*in[0]*in[2] + 72.0*in[0]*in[1]*in[2] - 27.0*in[0]*in[0] + 54.0*in[0]*in[0]*in[1] + 54.0*in[0]*in[0]*in[2] - 108.0*in[0]*in[0]*in[1]*in[2];
      out[19][1] = 0.0;
      out[19][2] = 0.0;

      out[20][0] = 0.0;
      out[20][1] = 9.0 - 18.0*in[0] - 36.0*in[1] - 18.0*in[2] + 72.0*in[0]*in[1] + 36.0*in[0]*in[2] + 72.0*in[1]*in[2] - 144.0*in[0]*in[1]*in[2] + 27.0*in[1]*in[1] - 54.0*in[1]*in[1]*in[0] - 54.0*in[1]*in[1]*in[2] + 108.0*in[0]*in[1]*in[1]*in[2];
      out[20][2] = 0.0;

      out[21][0] = 0.0;
      out[21][1] = -18.0*in[1] + 36.0*in[0]*in[1] + 36.0*in[1]*in[2] - 72.0*in[0]*in[1]*in[2] + 27.0*in[1]*in[1] - 54.0*in[0]*in[1]*in[1] - 54.0*in[1]*in[1]*in[2] + 108.0*in[0]*in[1]*in[1]*in[2];
      out[21][2] = 0.0;

      out[22][0] = 0.0;
      out[22][1] = 0.0;
      out[22][2] = 9.0 - 18.0*in[0] - 18.0*in[1] - 36.0*in[2] + 36.0*in[0]*in[1] + 72.0*in[0]*in[2] + 72.0*in[1]*in[2] - 144.0*in[0]*in[1]*in[2] + 27.0*in[2]*in[2] - 54.0*in[0]*in[2]*in[2] - 54.0*in[1]*in[2]*in[2] + 108.0*in[0]*in[1]*in[2]*in[2];

      out[23][0] = 0.0;
      out[23][1] = 0.0;
      out[23][2] = -18.0*in[2] + 36.0*in[0]*in[2] + 36.0*in[1]*in[2] - 72.0*in[0]*in[1]*in[2] + 27.0*in[2]*in[2] - 54.0*in[0]*in[2]*in[2] - 54.0*in[1]*in[2]*in[2] + 108.0*in[0]*in[1]*in[2]*in[2];

      out[24][0] = 96.0*in[0] - 144.0*in[0]*in[1] - 144.0*in[0]*in[2] + 216.0*in[0]*in[1]*in[2] - 96.0*in[0]*in[0] + 144.0*in[0]*in[0]*in[1] + 144.0*in[0]*in[0]*in[2] - 216.0*in[0]*in[0]*in[1]*in[2];
      out[24][1] = 0.0;
      out[24][2] = 0.0;

      out[25][0] = 0.0;
      out[25][1] = 96.0*in[1] - 144.0*in[0]*in[1] - 144.0*in[1]*in[2] + 216.0*in[0]*in[1]*in[2] - 96.0*in[1]*in[1] + 144.0*in[0]*in[1]*in[1] + 144.0*in[1]*in[1]*in[2] - 216.0*in[0]*in[1]*in[1]*in[2];
      out[25][2] = 0.0;

      out[26][0] = 0.0;
      out[26][1] = 0.0;
      out[26][2] = 96.0*in[2] - 144.0*in[0]*in[2] - 144.0*in[1]*in[2] + 216.0*in[0]*in[1]*in[2] - 96.0*in[2]*in[2] + 144.0*in[0]*in[2]*in[2] + 144.0*in[1]*in[2]*in[2] - 216.0*in[0]*in[1]*in[2]*in[2];

      out[27][0] = -144.0*in[0] + 288.0*in[0]*in[1] + 216.0*in[0]*in[2] - 432.0*in[0]*in[1]*in[2] + 144.0*in[0]*in[0] - 288.0*in[0]*in[0]*in[1] - 216.0*in[0]*in[0]*in[2] + 432.0*in[0]*in[0]*in[1]*in[2];
      out[27][1] = 0.0;
      out[27][2] = 0.0;

      out[28][0] = -144.0*in[0] + 216.0*in[0]*in[1] + 288.0*in[0]*in[2] - 432.0*in[0]*in[1]*in[2] + 144.0*in[0]*in[0] - 216.0*in[0]*in[0]*in[1] - 288.0*in[0]*in[0]*in[2] + 432.0*in[0]*in[0]*in[1]*in[2];
      out[28][1] = 0.0;
      out[28][2] = 0.0;

      out[29][0] = 0.0;
      out[29][1] = -144.0*in[1] + 288.0*in[0]*in[1] + 216.0*in[1]*in[2] - 432.0*in[0]*in[1]*in[2] + 144.0*in[1]*in[1] - 288.0*in[0]*in[1]*in[1] - 216.0*in[1]*in[1]*in[2] + 432.0*in[0]*in[1]*in[1]*in[2];
      out[29][2] = 0.0;

      out[30][0] = 0.0;
      out[30][1] = -144.0*in[1] + 216.0*in[0]*in[1] + 288.0*in[1]*in[2] - 432.0*in[0]*in[1]*in[2] + 144.0*in[1]*in[1] - 216.0*in[0]*in[1]*in[1] - 288.0*in[1]*in[1]*in[2] + 432.0*in[0]*in[1]*in[1]*in[2];
      out[30][2] = 0.0;

      out[31][0] = 0.0;
      out[31][1] = 0.0;
      out[31][2] =-144.0*in[2] + 288.0*in[0]*in[2] + 216.0*in[1]*in[2] - 432.0*in[0]*in[1]*in[2] + 144.0*in[2]*in[2] - 288.0*in[0]*in[2]*in[2] - 216.0*in[1]*in[2]*in[2] + 432.0*in[0]*in[1]*in[2]*in[2];

      out[32][0] = 0.0;
      out[32][1] = 0.0;
      out[32][2] = -144.0*in[2] + 216.0*in[0]*in[2] + 288.0*in[1]*in[2] - 432.0*in[0]*in[1]*in[2] + 144.0*in[2]*in[2] - 216.0*in[0]*in[2]*in[2] - 288.0*in[1]*in[2]*in[2] + 432.0*in[0]*in[1]*in[2]*in[2];

      out[33][0] = 216.0*in[0] - 432.0*in[0]*in[1] - 432.0*in[0]*in[2] + 864.0*in[0]*in[1]*in[2] - 216.0*in[0]*in[0] + 432.0*in[0]*in[0]*in[1] + 432.0*in[0]*in[0]*in[2] - 864.0*in[0]*in[0]*in[1]*in[2];
      out[33][1] = 0.0;
      out[33][2] = 0.0;

      out[34][0] = 0.0;
      out[34][1] = 216.0*in[1] - 432.0*in[0]*in[1] - 432.0*in[1]*in[2] + 864.0*in[0]*in[1]*in[2] - 216.0*in[1]*in[1] + 432.0*in[0]*in[1]*in[1] + 432.0*in[1]*in[1]*in[2] - 864.0*in[0]*in[1]*in[1]*in[2];
      out[34][2] = 0.0;

      out[35][0] = 0.0;
      out[35][1] = 0.0;
      out[35][2] = 216.0*in[2] - 432.0*in[0]*in[2] - 432.0*in[1]*in[2] + 864.0*in[0]*in[1]*in[2] - 216.0*in[2]*in[2] + 432.0*in[0]*in[2]*in[2] + 432.0*in[1]*in[2]*in[2] - 864.0*in[0]*in[1]*in[2]*in[2];
    }

    /**
     * \brief Evaluate Jacobian of all shape functions
     *
     * \param in Position
     * \param out return value
     */
    inline void evaluateJacobian (const typename Traits::DomainType& in,
                                  std::vector<typename Traits::JacobianType>& out) const
    {
      out.resize(36);

      out[0][0][0] = sign0*(-6.0*in[0] + 4);
      out[0][0][1] = 0;
      out[0][0][2] = 0;
      out[0][1][0] = 0;
      out[0][1][1] = 0;
      out[0][1][2] = 0;
      out[0][2][0] = 0;
      out[0][2][1] = 0;
      out[0][2][2] = 0;

      out[1][0][0] = sign1*(-2 + 6.0*in[0]);
      out[1][0][1] = 0;
      out[1][0][2] = 0;
      out[1][1][0] = 0;
      out[1][1][1] = 0;
      out[1][1][2] = 0;
      out[1][2][0] = 0;
      out[1][2][1] = 0;
      out[1][2][2] = 0;

      out[2][0][0] = 0;
      out[2][0][1] = 0;
      out[2][0][2] = 0;
      out[2][1][0] = 0;
      out[2][1][1] = sign2*(-6.0*in[1] + 4);
      out[2][1][2] = 0;
      out[2][2][0] = 0;
      out[2][2][1] = 0;
      out[2][2][2] = 0;

      out[3][0][0] = 0;
      out[3][0][1] = 0;
      out[3][0][2] = 0;
      out[3][1][0] = 0;
      out[3][1][1] = sign3*(6.0*in[1] - 2);
      out[3][1][2] = 0;
      out[3][2][0] = 0;
      out[3][2][1] = 0;
      out[3][2][2] = 0;

      out[4][0][0] = 0;
      out[4][0][1] = 0;
      out[4][0][2] = 0;
      out[4][1][0] = 0;
      out[4][1][1] = 0;
      out[4][1][2] = 0;
      out[4][2][0] = 0;
      out[4][2][1] = 0;
      out[4][2][2] = sign4*(-6.0*in[2] + 4);

      out[5][0][0] = 0;
      out[5][0][1] = 0;
      out[5][0][2] = 0;
      out[5][1][0] = 0;
      out[5][1][1] = 0;
      out[5][1][2] = 0;
      out[5][2][0] = 0;
      out[5][2][1] = 0;
      out[5][2][2] = sign5*(6.0*in[2] - 2);

      out[6][0][0] = -36.0*in[0]*in[1] + 18.0*in[0] + 24.0*in[1] - 12.0;
      out[6][0][1] = -18.0*in[0]*in[0] + 24.0*in[0] - 6;
      out[6][0][2] = 0.0;
      out[6][1][0] = 0.0;
      out[6][1][1] = 0.0;
      out[6][1][2] = 0.0;
      out[6][2][0] = 0.0;
      out[6][2][1] = 0.0;
      out[6][2][2] = 0.0;

      out[7][0][0]  = -36.0*in[0]*in[1] + 18.0*in[0] + 12.0*in[1] - 6.0;
      out[7][0][1] = -18.0*in[0]*in[0] + 12.0*in[0];
      out[7][0][2] = 0.0;
      out[7][1][0]  = 0.0;
      out[7][1][1] = 0.0;
      out[7][1][2] = 0.0;
      out[7][2][0]  = 0.0;
      out[7][2][1] = 0.0;
      out[7][2][2] = 0.0;

      out[8][0][0]  = 0.0;
      out[8][0][1]  = 0.0;
      out[8][0][2]  = 0.0;
      out[8][1][0]  = 18.0*in[1]*in[1] - 24.0*in[1] + 6.0;
      out[8][1][1]  = 36.0*in[0]*in[1] - 18.0*in[1] - 24.0*in[0] + 12.0;
      out[8][1][2]  = 0.0;
      out[8][2][0]  = 0.0;
      out[8][2][1]  = 0.0;
      out[8][2][2]  = 0.0;

      out[9][0][0]  = 0.0;
      out[9][0][1]  = 0.0;
      out[9][0][2]  = 0.0;
      out[9][1][0]  = -12.0*in[1] + 18.0*in[1]*in[1];
      out[9][1][1]  = 6.0 - 12.0*in[0] - 18.0*in[1] + 36.0*in[0]*in[1];;
      out[9][1][2]  = 0.0;
      out[9][2][0]  = 0.0;
      out[9][2][1]  = 0.0;
      out[9][2][2]  = 0.0;

      out[10][0][0] = 0.0;
      out[10][0][1] = 0.0;
      out[10][0][2] = 0.0;
      out[10][1][0] = 0.0;
      out[10][1][1] = 0.0;
      out[10][1][2] = 0.0;
      out[10][2][0] = 6.0 - 24.0*in[2] + 18.0*in[2]*in[2];
      out[10][2][1] = 0.0;
      out[10][2][2] = 12.0 - 24.0*in[0] - 18.0*in[2] + 36.0*in[0]*in[2];

      out[11][0][0] = 0.0;
      out[11][0][1] = 0.0;
      out[11][0][2] = 0.0;
      out[11][1][0] = 0.0;
      out[11][1][1] = 0.0;
      out[11][1][2] = 0.0;
      out[11][2][0] = -12.0*in[2] + 18.0*in[2]*in[2];
      out[11][2][1] = 0.0;
      out[11][2][2] = 6.0 - 12.0*in[0] - 18.0*in[2] + 36.0*in[0]*in[2];

      out[12][0][0] = -12.0 + 24.0*in[2] + 18.0*in[0] - 36.0*in[0]*in[2];
      out[12][0][1] = 0.0;
      out[12][0][2] = -6.0 + 24.0*in[0] - 18.0*in[0]*in[0];
      out[12][1][0] = 0.0;
      out[12][1][1] = 0.0;
      out[12][1][2] = 0.0;
      out[12][2][0] = 0.0;
      out[12][2][1] = 0.0;
      out[12][2][2] = 0.0;

      out[13][0][0] = -6.0 + 12.0*in[2] + 18.0*in[0] - 36.0*in[0]*in[2];
      out[13][0][1] = 0.0;
      out[13][0][2] = 12.0*in[0] - 18.0*in[0]*in[0];
      out[13][1][0] = 0.0;
      out[13][1][1] = 0.0;
      out[13][1][2] = 0.0;
      out[13][2][0] = 0.0;
      out[13][2][1] = 0.0;
      out[13][2][2] = 0.0;

      out[14][0][0] = 0.0;
      out[14][0][1] = 0.0;
      out[14][0][2] = 0.0;
      out[14][1][0] = 0.0;
      out[14][1][1] = -12.0 + 24.0*in[2] + 18.0*in[1] - 36.0*in[1]*in[2];
      out[14][1][2] = -6.0 + 24.0*in[1] - 18.0*in[1]*in[1];
      out[14][2][0] = 0.0;
      out[14][2][1] = 0.0;
      out[14][2][2] = 0.0;

      out[15][0][0] = 0.0;
      out[15][0][1] = 0.0;
      out[15][0][2] = 0.0;
      out[15][1][0] = 0.0;
      out[15][1][1] = -6.0 + 12.0*in[2] + 18.0*in[1] - 36.0*in[1]*in[2];
      out[15][1][2] = 12.0*in[1] - 18.0*in[1]*in[1];
      out[15][2][0] = 0.0;
      out[15][2][1] = 0.0;
      out[15][2][2] = 0.0;

      out[16][0][0] = 0.0;
      out[16][0][1] = 0.0;
      out[16][0][2] = 0.0;
      out[16][1][0] = 0.0;
      out[16][1][1] = 0.0;
      out[16][1][2] = 0.0;
      out[16][2][0] = 0.0;
      out[16][2][1] = 6.0 - 24.0*in[2] + 18.0*in[2]*in[2];
      out[16][2][2] = 12.0 - 24.0*in[1] - 18.0*in[2] + 36.0*in[1]*in[2];

      out[17][0][0] = 0.0;
      out[17][0][1] = 0.0;
      out[17][0][2] = 0.0;
      out[17][1][0] = 0.0;
      out[17][1][1] = 0.0;
      out[17][1][2] = 0.0;
      out[17][2][0] = 0.0;
      out[17][2][1] = -12.0*in[2] + 18.0*in[2]*in[2];
      out[17][2][2] =  6.0 - 12.0*in[1] - 18.0*in[2] + 36.0*in[1]*in[2];

      out[18][0][0] = 36.0 - 72.0*in[1] - 72.0*in[2] + 144.0*in[1]*in[2] - 54.0*in[0] + 108.0*in[0]*in[1] + 108.0*in[0]*in[2] - 216.0*in[0]*in[1]*in[2];
      out[18][0][1] = 18.0 - 72.0*in[0] - 36.0*in[2] + 144.0*in[0]*in[2] + 54.0*in[0]*in[0] - 108.0*in[0]*in[0]*in[2];
      out[18][0][2] = 18.0 - 72.0*in[0] - 36.0*in[1] + 144.0*in[0]*in[1] + 54.0*in[0]*in[0] - 108.0*in[0]*in[0]*in[1];
      out[18][1][0] = 0.0;
      out[18][1][1] = 0.0;
      out[18][1][2] = 0.0;
      out[18][2][0] = 0.0;
      out[18][2][1] = 0.0;
      out[18][2][2] = 0.0;

      out[19][0][0] = 18 - 36.0*in[1] - 36.0*in[2] + 72.0*in[1]*in[2] - 54.0*in[0] + 108.0*in[0]*in[1] + 108.0*in[0]*in[2] - 216.0*in[0]*in[1]*in[2];
      out[19][0][1] = -36.0*in[0] + 72.0*in[0]*in[2] + 54.0*in[0]*in[0] - 108.0*in[0]*in[0]*in[2];
      out[19][0][2] = -36.0*in[0] + 72.0*in[0]*in[1] + 54.0*in[0]*in[0] - 108.0*in[0]*in[0]*in[1];
      out[19][1][0] = 0.0;
      out[19][1][1] = 0.0;
      out[19][1][2] = 0.0;
      out[19][2][0] = 0.0;
      out[19][2][1] = 0.0;
      out[19][2][2] = 0.0;

      out[20][0][0] = 0.0;
      out[20][0][1] = 0.0;
      out[20][0][2] = 0.0;
      out[20][1][0] = -18.0 + 72.0*in[1] + 36.0*in[2] - 144.0*in[1]*in[2] - 54.0*in[1]*in[1] + 108.0*in[1]*in[1]*in[2];
      out[20][1][1] = -36.0 + 72.0*in[0] + 72.0*in[2] - 144.0*in[0]*in[2] + 54.0*in[1] - 108.0*in[1]*in[0] - 108.0*in[1]*in[2] + 216.0*in[0]*in[1]*in[2];
      out[20][1][2] = -18.0 + 36.0*in[0] + 72.0*in[1] - 144.0*in[0]*in[1] - 54.0*in[1]*in[1] + 108.0*in[0]*in[1]*in[1];
      out[20][2][0] = 0.0;
      out[20][2][1] = 0.0;
      out[20][2][2] = 0.0;

      out[21][0][0] = 0.0;
      out[21][0][1] = 0.0;
      out[21][0][2] = 0.0;
      out[21][1][0] = 36.0*in[1] - 72.0*in[1]*in[2] - 54.0*in[1]*in[1] + 108.0*in[1]*in[1]*in[2];
      out[21][1][1] = -18.0 + 36.0*in[0] + 36.0*in[2] - 72.0*in[0]*in[2] + 54.0*in[1] - 108.0*in[0]*in[1] - 108.0*in[1]*in[2] + 216.0*in[0]*in[1]*in[2];
      out[21][1][2] = 36.0*in[1] - 72.0*in[0]*in[1] - 54.0*in[1]*in[1] + 108.0*in[0]*in[1]*in[1];
      out[21][2][0] = 0.0;
      out[21][2][1] = 0.0;
      out[21][2][2] = 0.0;

      out[22][0][0] = 0.0;
      out[22][0][1] = 0.0;
      out[22][0][2] = 0.0;
      out[22][1][0] = 0.0;
      out[22][1][1] = 0.0;
      out[22][1][2] = 0.0;
      out[22][2][0] = -18.0 + 36.0*in[1] + 72.0*in[2] - 144.0*in[1]*in[2] - 54.0*in[2]*in[2] + 108.0*in[1]*in[2]*in[2];
      out[22][2][1] = -18.0 + 36.0*in[0] + 72.0*in[2] - 144.0*in[0]*in[2] - 54.0*in[2]*in[2] + 108.0*in[0]*in[2]*in[2];
      out[22][2][2] = -36.0 + 72.0*in[0] + 72.0*in[1] - 144.0*in[0]*in[1] + 54.0*in[2] - 108.0*in[0]*in[2] - 108.0*in[1]*in[2] + 216.0*in[0]*in[1]*in[2];

      out[23][0][0] = 0.0;
      out[23][0][1] = 0.0;
      out[23][0][2] = 0.0;
      out[23][1][0] = 0.0;
      out[23][1][1] = 0.0;
      out[23][1][2] = 0.0;
      out[23][2][0] = 36.0*in[2] - 72.0*in[1]*in[2] - 54.0*in[2]*in[2] + 108.0*in[1]*in[2]*in[2];
      out[23][2][1] = 36.0*in[2] - 72.0*in[0]*in[2] - 54.0*in[2]*in[2] + 108.0*in[0]*in[2]*in[2];
      out[23][2][2] = -18.0 + 36.0*in[0] + 36.0*in[1] - 72.0*in[0]*in[1] + 54.0*in[2] - 108.0*in[0]*in[2] - 108.0*in[1]*in[2] + 216.0*in[0]*in[1]*in[2];

      out[24][0][0] = 96.0 - 144.0*in[1] - 144.0*in[2] + 216.0*in[1]*in[2] - 192.0*in[0] + 288.0*in[0]*in[1] + 288.0*in[0]*in[2] - 432.0*in[0]*in[1]*in[2];
      out[24][0][1] = -144.0*in[0] + 216.0*in[0]*in[2] + 144.0*in[0]*in[0] - 216.0*in[0]*in[0]*in[2];
      out[24][0][2] = -144.0*in[0] + 216.0*in[0]*in[1] + 144.0*in[0]*in[0] - 216.0*in[0]*in[0]*in[1];
      out[24][1][0] = 0.0;
      out[24][1][1] = 0.0;
      out[24][1][2] = 0.0;
      out[24][2][0] = 0.0;
      out[24][2][1] = 0.0;
      out[24][2][2] = 0.0;

      out[25][0][0] = 0.0;
      out[25][0][1] = 0.0;
      out[25][0][2] = 0.0;
      out[25][1][0] = -144.0*in[1] + 216.0*in[1]*in[2] + 144.0*in[1]*in[1] - 216.0*in[1]*in[1]*in[2];
      out[25][1][1] = 96.0 - 144.0*in[0] - 144.0*in[2] + 216.0*in[0]*in[2] - 192.0*in[1] + 288.0*in[0]*in[1] + 288.0*in[1]*in[2] - 432.0*in[0]*in[1]*in[2];
      out[25][1][2] = -144.0*in[1] + 216.0*in[0]*in[1] + 144.0*in[1]*in[1] - 216.0*in[0]*in[1]*in[1];
      out[25][2][0] = 0.0;
      out[25][2][1] = 0.0;
      out[25][2][2] = 0.0;

      out[26][0][0] = 0.0;
      out[26][0][1] = 0.0;
      out[26][0][2] = 0.0;
      out[26][1][0] = 0.0;
      out[26][1][1] = 0.0;
      out[26][1][2] = 0.0;
      out[26][2][0] = -144.0*in[2] + 216.0*in[1]*in[2] + 144.0*in[2]*in[2] - 216.0*in[1]*in[2]*in[2];
      out[26][2][1] = -144.0*in[2] + 216.0*in[0]*in[2] + 144.0*in[2]*in[2] - 216.0*in[0]*in[2]*in[2];
      out[26][2][2] = 96.0 - 144.0*in[0] - 144.0*in[1] + 216.0*in[0]*in[1] - 192.0*in[2] + 288.0*in[0]*in[2] + 288.0*in[1]*in[2] - 432.0*in[0]*in[1]*in[2];

      out[27][0][0] = -144.0 + 288.0*in[1] + 216.0*in[2] - 432.0*in[1]*in[2] + 288.0*in[0] - 576.0*in[0]*in[1] - 432.0*in[0]*in[2] + 864.0*in[0]*in[1]*in[2];
      out[27][0][1] = 288.0*in[0] - 432.0*in[0]*in[2] - 288.0*in[0]*in[0] + 432.0*in[0]*in[0]*in[2];
      out[27][0][2] = 216.0*in[0] - 432.0*in[0]*in[1] - 216.0*in[0]*in[0] + 432.0*in[0]*in[0]*in[1];
      out[27][1][0] = 0.0;
      out[27][1][1] = 0.0;
      out[27][1][2] = 0.0;
      out[27][2][0] = 0.0;
      out[27][2][1] = 0.0;
      out[27][2][2] = 0.0;

      out[28][0][0] = -144.0 + 216.0*in[1] + 288.0*in[2] - 432.0*in[1]*in[2] + 288.0*in[0] - 432.0*in[0]*in[1] - 576.0*in[0]*in[2] + 864.0*in[0]*in[1]*in[2];
      out[28][0][1] = 216.0*in[0] - 432.0*in[0]*in[2] - 216.0*in[0]*in[0] + 432.0*in[0]*in[0]*in[2];
      out[28][0][2] = 288.0*in[0] - 432.0*in[0]*in[1] - 288.0*in[0]*in[0] + 432.0*in[0]*in[0]*in[1];
      out[28][1][0] = 0.0;
      out[28][1][1] = 0.0;
      out[28][1][2] = 0.0;
      out[28][2][0] = 0.0;
      out[28][2][1] = 0.0;
      out[28][2][2] = 0.0;

      out[29][0][0] = 0.0;
      out[29][0][1] = 0.0;
      out[29][0][2] = 0.0;
      out[29][1][0] = 288.0*in[1] - 432.0*in[1]*in[2] - 288.0*in[1]*in[1] + 432.0*in[1]*in[1]*in[2];
      out[29][1][1] = -144.0 + 288.0*in[0] + 216.0*in[2] - 432.0*in[0]*in[2] + 288.0*in[1] - 576.0*in[0]*in[1] - 432.0*in[1]*in[2] + 864.0*in[0]*in[1]*in[2];
      out[29][1][2] = 216.0*in[1] - 432.0*in[0]*in[1] - 216.0*in[1]*in[1] + 432.0*in[0]*in[1]*in[1];
      out[29][2][0] = 0.0;
      out[29][2][1] = 0.0;
      out[29][2][2] = 0.0;

      out[30][0][0] = 0.0;
      out[30][0][1] = 0.0;
      out[30][0][2] = 0.0;
      out[30][1][0] = 216.0*in[1] - 432.0*in[1]*in[2] - 216.0*in[1]*in[1] + 432.0*in[1]*in[1]*in[2];
      out[30][1][1] = -144.0 + 216.0*in[0] + 288.0*in[2] - 432.0*in[0]*in[2] + 288.0*in[1] - 432.0*in[0]*in[1] - 576.0*in[1]*in[2] + 864.0*in[0]*in[1]*in[2];
      out[30][1][2] = 288.0*in[1] - 432.0*in[0]*in[1] - 288.0*in[1]*in[1] + 432.0*in[0]*in[1]*in[1];
      out[30][2][0] = 0.0;
      out[30][2][1] = 0.0;
      out[30][2][2] = 0.0;

      out[31][0][0] = 0.0;
      out[31][0][1] = 0.0;
      out[31][0][2] = 0.0;
      out[31][1][0] = 0.0;
      out[31][1][1] = 0.0;
      out[31][1][2] = 0.0;
      out[31][2][0] = 288.0*in[2] - 432.0*in[1]*in[2] - 288.0*in[2]*in[2] + 432.0*in[1]*in[2]*in[2];
      out[31][2][1] = 216.0*in[2] - 432.0*in[0]*in[2] - 216.0*in[2]*in[2] + 432.0*in[0]*in[2]*in[2];
      out[31][2][2] = -144.0 + 288.0*in[0] + 216.0*in[1] - 432.0*in[0]*in[1] + 288.0*in[2] - 576.0*in[0]*in[2] - 432.0*in[1]*in[2] + 864.0*in[0]*in[1]*in[2];

      out[32][0][0] = 0.0;
      out[32][0][1] = 0.0;
      out[32][0][2] = 0.0;
      out[32][1][0] = 0.0;
      out[32][1][1] = 0.0;
      out[32][1][2] = 0.0;
      out[32][2][0] = 216.0*in[2] - 432.0*in[1]*in[2] - 216.0*in[2]*in[2] + 432.0*in[1]*in[2]*in[2];
      out[32][2][1] = 288.0*in[2] - 432.0*in[0]*in[2] - 288.0*in[2]*in[2] + 432.0*in[0]*in[2]*in[2];
      out[32][2][2] = -144.0 + 216.0*in[0] + 288.0*in[1] - 432.0*in[0]*in[1] + 288.0*in[2] - 432.0*in[0]*in[2] - 576.0*in[1]*in[2] + 864.0*in[0]*in[1]*in[2];

      out[33][0][0] = 216.0 - 432.0*in[1] - 432.0*in[2] + 864.0*in[1]*in[2] - 432.0*in[0] + 864.0*in[0]*in[1] + 864.0*in[0]*in[2] - 1728.0*in[0]*in[1]*in[2];
      out[33][0][1] = -432.0*in[0] + 864.0*in[0]*in[2] + 432.0*in[0]*in[0] - 864.0*in[0]*in[0]*in[2];
      out[33][0][2] = -432.0*in[0] + 864.0*in[0]*in[1] + 432.0*in[0]*in[0] - 864.0*in[0]*in[0]*in[1];
      out[33][1][0] = 0.0;
      out[33][1][1] = 0.0;
      out[33][1][2] = 0.0;
      out[33][2][0] = 0.0;
      out[33][2][1] = 0.0;
      out[33][2][2] = 0.0;

      out[34][0][0] = 0.0;
      out[34][0][1] = 0.0;
      out[34][0][2] = 0.0;
      out[34][1][0] = -432.0*in[1] + 864.0*in[1]*in[2] + 432.0*in[1]*in[1] - 864.0*in[1]*in[1]*in[2];
      out[34][1][1] = 216.0 - 432.0*in[0] - 432.0*in[2] + 864.0*in[0]*in[2] - 432.0*in[1] + 864.0*in[0]*in[1] + 864.0*in[1]*in[2] - 1728.0*in[0]*in[1]*in[2];
      out[34][1][2] = -432.0*in[1] + 864.0*in[0]*in[1] + 432.0*in[1]*in[1] - 864.0*in[0]*in[1]*in[1];
      out[34][2][0] = 0.0;
      out[34][2][1] = 0.0;
      out[34][2][2] = 0.0;

      out[35][0][0] = 0.0;
      out[35][0][1] = 0.0;
      out[35][0][2] = 0.0;
      out[35][1][0] = 0.0;
      out[35][1][1] = 0.0;
      out[35][1][2] = 0.0;
      out[35][2][0] = -432.0*in[2] + 864.0*in[1]*in[2] + 432.0*in[2]*in[2] - 864.0*in[1]*in[2]*in[2];
      out[35][2][1] = -432.0*in[2] + 864.0*in[0]*in[2] + 432.0*in[2]*in[2] - 864.0*in[0]*in[2]*in[2];
      out[35][2][2] = 216.0 - 432.0*in[0] - 432.0*in[1] + 864.0*in[0]*in[1] - 432.0*in[2] + 864.0*in[0]*in[2] + 864.0*in[1]*in[2] - 1728.0*in[0]*in[1]*in[2];
    }

    //! \brief Evaluate partial derivatives of all shape functions
    void partial (const std::array<unsigned int, 3>& order,
                  const typename Traits::DomainType& in,         // position
                  std::vector<typename Traits::RangeType>& out) const      // return value
    {
      auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
      if (totalOrder == 0) {
        evaluateFunction(in, out);
      } else {
        DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
      }
    }

    //! \brief Polynomial order of the shape functions
    unsigned int order () const
    {
      return 3;
    }

  private:
    R sign0, sign1, sign2, sign3, sign4, sign5;
  };
}
#endif // DUNE_LOCALFUNCTIONS_RAVIARTTHOMAS1_CUBE3D_LOCALBASIS_HH
